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ABSTRACT 

We present three-dimensional, nonrelativistic, hydrodynamic simulations of bow 
shocks in pulsar wind nebulae. The simulations are performed for a range of initial 
and boundary conditions to quantify the degree of asymmetry produced by latitudinal 
variations in the momentum flux of the pulsar wind, radiative cooling in the postshock 
flow, and density gradients in the interstellar medium (ISM). We find that the bow 
shock is stable even when travelling through a strong ISM gradient. We demonstrate 
how the shape of the bow shock changes when the pulsar encounters density variations 
in the ISM. We show that a density wall can account for the peculiar bow shock shapes 
of the nebulae around PSR J2124-3358 and PSR B0740-28. A wall produces kinks in 
the shock, whereas a smooth ISM density gradient tilts the shock. We conclude that 
the anisotropy of the wind momentum flux alone cannot explain the observed bow 
shock morphologies but it is instead necessary to take into account external effects. 
We show that the analytic (single layer, thin shell) solution is a good approximation 
when the momentum flux is anisotropic, fails for a steep ISM density gradient, and ap- 
proaches the numerical solution for efficient cooling. We provide analytic expressions 
for the latitudinal dependence of a vacuum-dipole wind and the associated shock 
shape, and compare the results to a split-monopole wind. We find that we are unable 
to distinguish between these two wind models purely from the bow shock morphology. 
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1 INTRODUCTION 

After a pulsar escapes its supernova remnant, it typi- 
cally travels supersonically through the interstellar medium 
(ISM), sometimes with speeds in excess of 1000 km s _1 
Ichatteriee et aljl2005h . At the same time, pulsars emit 
highly relativistic winds, probably consisting of electron- 
p osit ron pairs, with bulk Lorentz factor between 10 2 an d 
10 6 jKennel fc Coronitil Il984al : ISpitkovskv fc Aronsl l2004) . 
The interaction of the pulsar wind with t he ISM produces a 
characteristic multilayer sh ock structure feucciantinil2002l 
van der Swaluw et alj|2003l) . 

The multilayer structure of pulsar wind nebulae (PWN) 
gives rise to a variety of observa b le feat ures at different wave- 
lengths; see iGaensler fc Slanel i2006f) for a recent review. 
The outer bow shock excites neutral hydrogen atoms either 
collisionally or by charge exchang e, which then de-excite 
and emit optical Ha radiation dBucciantini fc Bandieral 
l200lT> . The microphysics of these excitation processes 



was treated in the context of s upernova remn a nts b y 
IChevalier fc Raymond! lll978l) and IChevalier et alJ <198Cft . 
while iGhavamiaii^t^l ll200lh performed detailed calcula- 
tions of the ionization structure and optical spectra pro- 
duced by non-radiative supernova remnants in partially neu- 
tral gas. Presently, six PWN have been detected in the opti- 
cal band iGaensler fc Slanel2006l : lKaspi etlii 12004). namely 
aroun d the pulsars PS R B 1957+20 faulkarni fc Hester! 
Il988ft . PSR B2224+65 (ICordes et ail Jl993l). also calle d 
the Guitar Nebula PSR B0740-28 i Jones et alJ 120021) . 
PSR J0437-4715 feell et alJ Il995l) . RX J1856. 5-3754 
I van Kerkwiik fc Kulkarnil l200ll) .~and PSR J2124-3358 
i Gaensler et al.ll2002l) ~ 



Pulsar bow shocks directly probe the conditions in the 
local ISM. T he size of the bow shock scales with the s tand- 
off distance JChatteriee fc Corded liool IWilkinlll996l) . de- 
fined as the point where the ISM ram pressure balances 
the wind ram pressure. As the wind ram pressure is given 
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directly by the (measurable) sp in-down luminosity of the 
pulsar for both young pulsars iKermel_j£_Co£orritj| h 
and recycled millisecond pulsars (|Stappers et alJl2003h . the 
size of the bow shock can put an upper limit on the ISM 
density, given a measurement of the pulsar's proper mo- 
tion. Furthermore, the morphology of the shock gives de- 
tailed information about the inhomogeneity of the ISM. 
Density gradients have been proposed to explain the peculiar 
shapes of the bow s hocks around the pulsars PSR B0740— 28 
] Jones et alJ l2002h and PSR J2124-3358 llGaensler et alJ 
l2002l) . IChatteriee et ail fcOOfJ) recently showed that the kink 
in the PWN around PSR J2124-3358 can be explained if 
the pulsar is travelling through a density discontinuity in 
the ISM (a "wall"). 

Our theoretical understanding of the wind acceleration 
mechanism, in particular the an gular distribution of i ts mo- 
mentum flux, is still incomplete. iGoldreich fc Julian! lll969T) 
showed that a pulsar has an extended magnetosphere filled 
with charged p articles that ar e accelerated into a surround- 
ing wave zone <Coronitill99ch . Although the nature of the 
energy transport in the pulsar wind is still an unsolved ques- 
tio n probably governed by no n-ideal magnetohydrodynam- 
K'S llMelatos fc Melrose! Il996l) . all models agree that it is 
dominated by the Poynting flux near the star. However, ob - 
servations of the Cr ab nebula |K ennel fc CoronitH Il984al) 
and PSR B1509-58 iGaensler et al]l2002h imply a value for 
the magnetisation parameter, a, defined as the ratio of the 
Poynting flux to the kinetic energy flux, satisfying a < 0.01 
at the wind termination shock. In other words, a transi- 
tion from high a to low a must occur between the magne- 
tosphere and the termination shock, a problem commonly 
referred to as the a paradox. While it is generally assumed 
that the angular dependenc e of the momentum flu x is pre- 
served during this transition Zanna et al I2004I) . the ex- 
act dependence is unknown. T he split-monopole model pre- 
dicts a mainly equatorial flux feogovaloylll999l: I Bucciantinil 
2006 ) while wave-like dipole models iM^iatos^^Melrose 



1996) predict a more complicated distribution. The study 
of the termination shock structure, including it s variability 
iSpitkovskv fc Aronsl[2004l : iMelatos et al] 120051) . thus pro- 
vides valuable information about the angular distribution of 
the wind momentum flux and hence the electrodynamics of 
the pulsar magnetosphere. 

PWN also e mit synchrotron radiation at radio and X - 
ray wavelengths llGaensler fc Slanell2006l : iKaspi et al]l2004Tl 
from wind particles which are accelerated in the inner 
shock and subsequently gyrate in the nebular magneti c 
field (iKennel fc Coronitill984bl: ISpitkovskv fc Aronslll999l) . 
If the pulsar moves supersonically, the (termin ation) shock 
is strongly confined by the ISM r am pressure feucciantinil 
120021 : Ivan der Swaluw et alj|2003h and elongates opposite 
the direction of motion. X-ray-emitting particles with short 
synchrotron lifetimes directly probe the inner flow structure 
of the PWN. With the help of com bined X-ray and radio 
observations, iGaensler et al] i2004f) identified the termina- 
tion shock and post shock flow in the PWN around PSR 
J1747-2958, the Mouse Nebula. They found that the data 
are consistent with an isotropic momentum flux distribu- 
tion, a statemen t that holds for some other bow shocks, e.g. 
PSR B 2224-65 jChatteriee fc Cordes^Ooll . but not all, e.g. 
IC433 jGaensler fc Slaneil2006ft ~ 

Before we can use PWN to probe the ISM, a precise 



knowledge of the inner flow structure and the respons e of the 
system to external influences is important. Although lWilkinl 
lll996l) found an analytic solution for the shape of the bow 
shock in the limit of a thin shock, the complexity of the prob- 
lem g enerally requires a numerical treatment. iBucciantinil 
i2002fl performed hydrodynamic, 2.5-dimensional, cylindri- 
cally symmetric simulations to explore the validity of a 
two-shell approximation to the pr oblem and elucidated 
the in ner, multilayer shock structure. Ivan der Swaluw et al] 
]2003h modelled a pulsar passing through its own super- 
nova remnant and found that the bow shock nebula remains 
undisrupted. 

Of course, the validity of a hydrodynamic approach is 
limited. To compute synchrotron maps, magnetic fields must 
be included. R elativistic magnetohydrodynami c (RM HD) 
simulations by iKomissarov fc Lvubarskvl ll2003l . 12004]) ex- 
plain the peculiar jet-like emission in the Crab nebula by 
magnetic collimation of the back fl ow of an anisotropic pu l- 
sar wind, a result confirmed later bv lDel Zanna eta]] (12004). 
RMH D simulations of bow sho cks around fast moving pul - 
sars bv lBucciantini et al] ]2005l) and lBogovalov et al.l (120051) 
assume cylindrical symmetry and an isotropic pulsar wind. 

The chief new contribution of this paper is to relax the 
assumption of cylindrical symmetry. We perform fully three- 
dimensional, hydrodynamic simulations of PWN, focusing 
on observable features which can be compared directly to 
observations. Although pulsar winds are generally expected 
to be symmetric around the pulsar rotation axis (even if 
the momentum flux varies strongly with latitude), there are 
external factors, like an ISM density gradient, which de- 
stroy the cylindrical symmetry. The existence of two very 
different speeds, namely that of the relativistic wind and 
the pulsar proper motion, presents a significant numerical 
challenge in three dimensions, which we overcome by us- 
ing a parallelised hydrodynamic code, with adaptive-mesh 
Godunov-type shock-capturing capabilities. 

After describing the numerical method in section [5] wc 
present the results of our simulations for different types of 
asymmetry, namely anisotropic momentum flux (section 0, 
the effects of uniform cooling (section , an ISM density 
gradient (section |5J, and an ISM barrier or "wall" (section 
|SJ. We conclude by discussing what information can be de- 
duced about the ISM and the pulsar wind when the results 
are combined with observations in section |7| 



2 NUMERICAL METHOD 
2.1 Hydrodynamic model 

The simulations are performed using Flash ]Frvxell et al] 
I2000T) . a parallelized hydrodynamic solver based on the 
second-order piecewise parabolic method (PPM). Flash 
solves the inviscid hydrodynamic equations in conservative 
form, 



dp 
dt 

d(pv 



dt 
8(pE) 



dt 



+ V ■ (pv) , 

+ V ■ (pvv) + VP, 
+ V ■ [{pE + P) v] 
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together with an ideal gas equation of state, 
P = (j-l)pe, 



(4) 



where p, P, e and v denote the density, pressure, inter- 
nal energy per unit mass, and velocity, respectively, and 
E = e + | -tf | 2 /2 is the total energy per unit mass. We take 
advantage of the multifluid capability of Flash to treat the 
pulsar wind and the interstellar medium as fluids with differ- 
ent adiabatic indices 71 and 72, respectively. The weighted 
average adiabatic index 7 can then be computed from 

1 X 1 , X 2 



7 " 



7i 



72 



where Xi is the mass fraction of each fluid, advected accord- 
ing to 



+ V ■ (pXiv) 



0. 



(6) 



d(pXj) 
dt 

We do not consider reactive flows, so there is only a small 
amount of numerical mixing between the two fluids. For the 
problem studied here, we find empirically that the results 
are qualitatively the same for 71 = 72 and 71 7^ 72 (see 
section |3~^1 and therefore take 71 =72 = 5/3 in most of our 
simulations, using the multifluid capability only to trace the 
contact discontinuity. 

Fl ash manages the a daptive mesh with the Paramesh 
library dOlson et al.lll999t) . The mesh is refined and coars- 
ened in response to the secon d-order error i n the dynam- 
ical variables (for details, see lLoehnerlll98^) . We perform 
our simulations on a three-dimensional Cartesian grid with 
8x8x8 cells per block and a maximum refinement level 
of five nested blocks, giving a maximum resolution of 128 3 
cells. 

We are limited by processing capacity to simulations 
lasting ~ 15ro, where to is the time for the ISM to cross the 
integration volume. Unfortunately, this is not always suffi- 
cient to reach a genuine steady state. We occasionally en- 
counter ripples in the shock [e.g. Figure (right), discussed 
further in section [3] which are numerical artifacts which di- 
minish as the simulation progresses. 



2.2 Boundary and initial conditions 

The simulations are performed in the pulsar rest frame. Fig- 
ure Q depicts the pulsar and the orientation of its rotation 
axis CI with respect to the velocity vector v m of the ISM 
in this frame. The momentum flux of the pulsar wind is as- 
sumed to be cylindrically symmetric about tt, in keeping 
with simplified t h eoretical models like the split m onopolc 
feogovalovl Il999t iKomissarov fc Lvubarskvl 120031) or vac- 
uum dipole llMelatos fc Mekosdll996UMelatodll997ll2002lh 
In the pulsar's rest frame, the ISM appears as a steady 
wind that enters the integration volume from the upper 
boundary z = 2 ma x, with velocity v m = — v p = —4 x I0 7 
cm s -1 z, where v H is the pulsar's velocity in the observer's 



frame (iHobbs et al .120051) . density 



p m = (0.33-1.30) x 10" 23 
g cm -3 , and specific internal energy e m = 1.18 x 10 12 erg 
g _1 . The sound speed c m = 1.15 x 10 6 cm s _1 and Mach 
number M m — v m /c m — 35 are typical of the parameters 
inferred for a variety o f observed PWN llKaspi et al.ll20o"il ; 
iBucciantini et al. 2005). All boundaries except z = z ma x act 
as outflow boundaries, where the values of all the dynamical 




Figure 1. Orientation of the pulsar, wall, and Cartesian coor- 
dinate axes in the pulsar's rest frame. The pulsar is located at 
(0, 0, z p ) and the interstellar medium moves parallel to the z axis 
with a speed v m . The symmetry axis f2 of the wind (assumed to 
be the pulsar's rotation axis) is defined by its polar angle A and 
azimuthal angle <j>. The wall lies perpendicular to the y-z plane, 
making an angle a with respect to the z axis. 



variables are equalized in the boundary cells and integration 
region (zero gradient condition). These boundary conditions 
are justified provided that the flow speed across the bound- 
aries exceeds the sound speed thus preventing a casual con- 
nection of the regions inside and outside the computational 
domain. We check this assumption a posteriori in section 
13.21 Initially, we set p — p m , v = (0, 0, — v m ), and e = e m in 
every cell. 

We allow for the possibility that the pulsar travels into 
a density gradient. This may be either a smooth gradient 
perpendicular to the pulsar's direction of motion (section 
|SJ or an extended ridge of high density material making an 
angle a with v m , which we call a "wall", also depicted in 
figureQ(section|5|l. The boundary conditions corresponding 
to these scenarios are defined in section |S] and HJ 

The simulation parameters for the different models are 
listed in tabled 



3 ANISOTROPIC WIND 

Although the electrodynamics of pulsar wind acceleration 
and collimation rem a ins unsolved, observa t ions of PWN 
jHelfand et alj l200ll: iGaensler etail 120021 iRoberts et all 
l2005fl as well as theoretical studies of split-mono pole 
jBogovakwlll999l: IKomissarov fc Lvubarskvl2'ool 120041) and 
wave- like dipole jMelatos fc Melrosel Il99& iMelatosI Il997i 
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2002) outflows seem to favor anisotropic momentum and 
energy flux distributions, where the wind is cylindrically 
symmetric about fl but varies in strength as a function 
of latitude. Letting 6 denote the colatitude measured rel- 
ative to f2, the various models predict the following mo- 
mentum flux di stributions : g w (8) = sin 2 6 + 1/cr for the 
split monopole l|Bogovalovlll999t iKomissarov fc Lvubarskvl 
120031: lAronsll2004h . a,„(0) = cos 2 6 + 1 for the point dipole, 
and g w (0) = 2 cos 4 6 + 6 cos 2 8 + 5 for the extended dipole 
see also appendix [XJ. The normalized mo- 
mentum flux g m {9) is defined precisely in appendix 151 



3.1 Numerical implementation 

We implement the wind anisotropy by varying the wind den- 
sity p w according to 



Pw(0) = po, w g w {Q) = po,w (co + c 2 < 



(7) 



where co, C2 and p are constants, while leaving the wind 
velocity v w constant. Strictly speaking, this is unrealistic, 
because it is likely that v„ varies with 9 too in a true pulsar 
wind. However, it is expected to b e an excellent approxi- 
mation because Iwilkinl jl996l . I2000F) showed that the shape 
of the bow shock is a function of the momentum flux p w v^ 
only, not p m and v w separately, at least in the thin-shell 
limit of a radiative shock. We deliberately choose equation 
to h ave exactly the same form as in the analytic theory 
(Wilkin 1996, 2000) to assist comparison below. 

In equation 0, Co and C2 parametrize the anisotropy 
of the wind, while p can be used to model jet-like outflows 
(e.g. p ^ 4). For p = 2, eq uation is equivalent to equa- 
tion (109) of lWilkinl feOOOh The split-monopole model cor- 
respo nds to p = 2 and C2 = —1 iKomissarov fc Lvubarskvl 
120031) . Note that normalisation requires Co = 1 — C2/3 for 
p = 2 and co = 1 — C2/5 for p = 4. The analytic solutions 
for these two cases are presented in appendix 1X1 

Although the pulsar wind is ultrarelativistic 
[e.g. Lorentz fact o r 10 2 — 10 6 in the Crab Nebul a; 
iKennel fc Coronitil dl984aTl ; ISpitkovskv fc Aronsl <2004) 1. 
the shape is determined by th e ratio of the momentum 
fluxes alone dWilkinlll996l . 120001) . We can therefore expect 
our model to accurately describe the flow structure, as long 
as the value of the mean flux po.wV^ = 4.6 x 10~ 8 g cm -1 
s~ 2 is realistic (which it is), even thou gh we have v w <S c . 
Thi s approach was also adopt e d bv IBu cciantinil ( 2002 ) 
and Ivan der Swaluw et al.1 (12003ft . iBucciantini et all <2005li 



performed RMHD simulations of PWN and a comparison 
of our results with theirs does not reveal any specifically 
relativistic effects; RMHD simulations with low a agree 
qualitatively with non-relativistic hydrodynamical simula- 
tions. In other words, as long as we choose p w and v w such 
that pwV^j = E/iirr w c, where E is the pulsar's spin-down 
luminosity (and r w is defined in the next paragraph), we 
get the same final result. 1 



1 It is important to bear in mind that this is not strictly true 
in an adiabatic (thick-shell) shock iLuo et alj|l99(jl or when the 
two winds contain different particle species (see Appendix B of 
iMelatos et al]|l995ft . where p w and v m control the shape of the 
shock separately and one should use the relativistic expression 
n w f w mc 2 v w for the momentum flux, where n w is the number 



The pulsar wind has density po, w = 4.6 x 10 -24 g cm~ 3 , 
specific internal energy t w = 5.9 x 10 12 erg g _1 , and velocity 
v w = v w e r , where e r is the radial unit vector and v w — 10 8 
cm s _1 . This particular choice gives a sound speed c w = 
2.5 x 10 6 cm s _1 and Mach number M w = 40. The pulsar 
wind is launched by restoring the dynamical variables p w , v„ 
and e w in a spherical region with radius r m = 6.2 x 10 15 cm 
around r p = (0,0, z p ) (figure to their initial values after 
every timestep. 



3.2 Multilayer shock structure 

Figure [5] shows the y-z section of model A, a typical 
anisotropic wind with p = 2 and C2 = 3 in which fl is 
tilted by 45°with respect to v p . This wind is pole-dominated, 
with zero momentum flux at the equator, contrary to what 
is inferred from observation jGaensler fc Sland l2006ft . We 
postpone a detailed comparison between a pole-dominated 
and an equatorially-dominated outflow to section 13.41 and 
discuss here only the shock structure, which is common to 
both cases. 

All our simulations show the ch aracteristic multilayer 
struct ure in figureEl des c ribed also bv lvan der Swaluw et alJ 
teOOSl) and IBucciantinil i2002f) . The pulsar wind inflates a 
wind cavity that is enclosed by a termination shock (TS). 
The location of the TS behind the pu lsar is determined 
by the thermal pressure i n the ISM iBucciantinil 1200 2l ; 
van der Swaluw et al J 1200 A . Ahead of the pulsar, the ter- 
mination shock is confined by the ram pressure of the ISM. 
The wind material is thermalised at the TS and fills a cylin- 
drical postshock region (PS), separated from the shocked 
ISM by a contact discontinuity (CD). 

In the case of a thi n shock ; the s hape of the CD can be 
described analytically JWilkinll200(i see also appendix lAt. 
with a global length scale Ro which is computed by equating 
the ram pressures of the ISM and the wind: 



Ro 



E 



4irp m VmV w 



1/2 



(8) 



The analytic solution is depicted by a solid curve (scaled to 
match the CD) and by a dotted curve (scaled to match the 
BS) in figurc|21 As discussed below (section EjSt i it describes 
the location of the CD and the BS, respectively, reasonably 
well. 

For an isotropic wind, _Ro equals the standoff distance, 
defined as the distance from the pulsar to the intersection 
point of the pulsar's velocity vector and the bow shock apex. 
For an anisotropic wind, however, the standoff distance also 
depends on the exact distribution of moment um flux and R p 
is merely an overall length-scale parameter iWilkinll2000l) . 

We can observe Kelvin-Helmholtz (KH) instabilites at 
the CD, most prominently in the lower-right corner of fig- 
ure|2 In the absence of gravity, KH instabilities occur on all 
wavelengths, down to the grid resolution, with growth rate 
merely proportional to wavelength squared. Sure enough, 
the length-scale of the ripples in figure|5|is indeed compara- 
ble to the grid resolution (A^ « 0.05i? , Ay ~ 0.125i? () ). 

Figure|3|compares the results for model A with 71 =5/3 



density, 7^ is the Lorentz factor, and m is the particle mass, as 
well as a relativistic code. 
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Table 1. Model parameters. The wind momentum flux is proportional to 1 — C2/3 + ci cos p 9, where 6 is the colatitude measured with 
respect to the pulsar spin axis f2, and A and <j> define the orientation of this axis (figure and section ^0 is the standoff distance as 
measured from the simulation output, / stands for the cooling efficiency (section 0, and H is the scale-height of the ambient density 
gradient (section |SJ. The wall (section [SJ is described by the density contrast n = p m /pm, the angle a, and the distance over which the 
density rises (A r ) and declines (A d ). The dots indicate parameters that are not used in that simulation. 
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Figure 4. The ratio of total velocity v = (v 2 , + V 2 + u|) lj ' 2 to 
adiabatic sound speed c s = ("/p/p) 1 ' 2 in a y-2 section of model A. 
The colors denote v/c s in logarithmic scale, as indicated by the 
color bar. The bulk flow is highly supersonic along, and nearly 
everywhere within, the boundaries. The solid lines denote the 
contour levels where v/c s = 1. 



(left) and 71 = 4/3 (right), while 72 = 5/3 in both cases. As 
expected, the overall shock morphology is weakly affected by 
adopting the adiabatic index 71 = 4/3 for a relativistic gas, 
except for the ripples which arise because the simulation has 
not yet reached a steady state. 

Figure 2] displays the Mach number v / c s [where c s = 
(iP/p) 1 ^ 2 is the adiabatic sound speed] as a function of 
position in a y-z section of model A. The flow is super- 
sonic nearly everywhere along the boundaries and on the 
pulsar side of the BS. The only exception is a few small 
regions of subsonic flow in the PS (solid contours in fig- 
ure and we do not expect the back reaction from these 
subsonic inclusions to affect the flow structure. In order to 
verify this claim, we perform one simulation with the pa- 



rameters of model A and a larger computational domain 

(•Cmax 3-miii — |/max Ulain — 2 ma x ^min — 21.3-RfJ, Com- 
pared to ^max ^min — |/max — 2 m ax ^min — 10.6-Zirj 

in figure |5J and compare the results in figure [S] We restrict 
the runtime to t = 3.38ro (where to is the time for the 
ISM to cross the visible domain) to keep the computation 
tractable. Although there are slight changes in the shock 
structure, which can be attributed to the different resolu- 
tion of the wind region, we do not see any evidence for back 
reaction (e.g., matter or sound waves travelling into the vis- 
ible domain). 

3.3 Shock asymmetries 

The analytic solution 1A21 (solid curve in figure approxi- 
mates the shape of the CD (dashed curve) reasonably well. 
Here, and in all the following plots, the CD in the simu- 
lation is defined by the condition X2 — 0.9 [see equation 
0]. Plainly, the CD in figure |2] is asymmetric. The kink 
that appears in the analytic solution at the latitude where 
PwV^ is a minimum (6 = tt/2, perpendicular to f2 in the fig- 
ure) is less distinct in the simulations, because the thermal 
pressure in the wind and the ISM smoothes out sharp fea- 
tures. Hence, even a highly anisotropic wind does not lead to 
large, easily observable, kink-like asymmetries in the CD, an 
important (if slightly disappointing) point to bear in mind 
when interpreting PWN observational data. 

The point of the BS closest to the pulsar is no longer 
along v p but depends also on the latitude where the mo- 
mentum flux is a minimum. However, the BS itself is also 
smoothed by thermal pressure so it is not a good gauge of 
the underlying anisotropy. Equally, it is difficult to infer the 
density of the ambient ISM from the location of the stagna- 
tion point, defined as the intersection point of v p with the 
BS, since the orientation of ft cannot be inferred uniquely 
from the shape of the BS. 
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Figure 2. section of model A. The colors denote the density in logarithmic scale, as indicated by the color bar. The characteristic 
multilayer structure, including the bow shock (BS), termination shock (TS), contact discontinuity (CD), and postshock region (PS), 
can be seen clearly. The solid curve shows the analytic approximation IA2I . with Rq scaled to match the CD (dashed curve), measured 
empirically by tracing the boundary between the two fluids in Flash. The dotted curve is the analytic solution scaled to match the BS. 
The projected symmetry axis and the wind momentum flux, drawn as a polar plot, appear in the top left corner of the figure. 



The asymmetry of the BS in figure[3] which is more ex- 
tended towards the right side of the figure, is c aused by the 
excess momentum flux along — v m on this side. iBucciantina 
i2002l) showed that the analytic formula 1AH accurately de- 
scribes the CD for isotropic momentum flux. We confirm 
here that the analogous formula 1A2II accurately describes 
the CD for anisotropic momentum flux. But, as the wind 
ram pressure decreases downstream from the apex, the ther- 
mal pressure becomes increasingly important and the ana- 
lytic solution breaks down. Nevertheless, the asymmetry is 
preserved. Likewise, the analytic solution for the BS (dotted 
curve) deviates more strongly as one moves further down- 
stream from the apex. Again, this result is expected, since 
the analytic approximation relies on mechanical momentum 
flux balance and breaks down as the thermal pressure be- 
comes increasingly important. 

Momentum flux anisotropy affects the BS weakly, but 
it affects the structure of the inner flow strongly. A striking 
feature of figure [5] is the shape of the TS. It is highly asym- 
metric and elongated, roughly aligned with tt. As the rear 
surface of the TS lies well inside the BS, its position is de- 
termined by the balance between the wind ram pressure and 
the thermal pressure in the PS region, which is found em- 



pirically to be roughly the same as in the ISM feucciantinil 
120021) . This, combined with the small momentum flux at 
= — 7r/2 (running diagonally from the upper left to bot- 
tom right in figure |5J, explains why the side surfaces of the 
TS are confined tightly, and hence its characteristic convex 
shape. 

We now contrast the above results with a system that 
is axisymmetric around z. Figure |S| shows the y-z section of 
model B, which is the same as model A, except that ft is 
now perpendicular to v m . Not surprisingly, the BS and CD 
are axisymmetric. The analytic formula HA2t describes the 
CD reasonably well, although the BS deviates slightly near 
the apex for the pressure-related reasons discussed above. 
The butterfly-like shape of the TS directly reflects the mo- 
mentum flux and is a good observational probe in this sit- 
uation. Note that the CD and BS nearly touch the pulsar 
at — —n/2 (where the momentum flux is zero). This illus- 
trates that an estimate of p m based on equation JSJ and the 
assumption of an isotropic wind is too low. Note again the 
KH instabilites at the CD and the blobs of ISM matter in 
the PS region, which are numerical artifacts. 

An example of the difficulties arising f rom the degen- 
eracy of the problem is PSR J2124— 3358 iGaensler et alJ 
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Figure 3. y-z section of model A for 71 = 5/3 (left) and 71 = 4/3 (right), while 72 = 5/3. The shock morphology is the same. The 
system has not yet reached steady state, as can be seen by the ripples in the BS and CD. 
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Figure 5. y-z sections of model A with a computational domain of a; ma x — a; m in = ?/max — J/min = Zmax — z m in = 10.6-Ro (left) and the 
same model with a domain of x max — aJmin = Jmax — y m in = Zmax — z m in = 21.3-Ro (right, zoomed-in to allow comparison with the left 
figure) at t = 3.38ro and t = 3.35to, respectively. The colors denote the density in logarithmic scale as in figurc[5] The shock morphology 
shows slight differences that can be attributed to the different evolution of the KH instabilities (due to the different resolution of the 
wind region and the time when the snapshot is taken). We do not observe any back reaction effects in the larger computational domain. 



2002). The PWN around this pulsar sh ows a clear asymme- 
try near the apex. lOaensrer et alJ l|2002F) showed that several 
combinations of an ISM density gradient (perpendicular and 
parallel to v m ), an ISM bulk flow, and an anisotropic pul- 
sar wind can explain the observed anisotropy, but only in 
combination, not alone. 

Figure |7| demonstrates what the BS looks like when the 



wind is highly collimated along ft, i.e. a jet (g w oc cos 4 #). 
Basically, the situation is similar to model A (figure |5J: the 
analytic formula correctly predicts the shape of the CD. 
However, the kink at the latitude of the momentum flux 
minimum is less visible than in figures|5|and|S| Moreover, the 
wind cavity is smaller than in figure [5] The TS approaches 
the pulsar more closely from behind, where there is less wind 




Figure 6. y-z section of model B. The colors denote the density in 
logarithmic scale, as defined in figure^] The solid curve shows the 
theoretical solution IA2I . scaled to match the CD (dashed curve). 
The solid and dashed curves overlap almost perfectly. The symbol 
in the top left indicates the projected wind symmetry axis and 
momentum flux (polar plot), as in the caption of figure |3 Note 
the butterfly-like TS and the BS almost touching the pulsar. 
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Figure 7. y-z section of model C. The colors denote the density in 
logarithmic scale, as defined in figure^! The symbol in the top left 
indicates the projected wind symmetry axis and momentum flux 
(polar plot), as in the caption of figure The solid curve shows 
the theoretical solution IA11I . scaled to match the CD (dashed 
curve). The momentum flux in the wind is highly collimated, with 
g w oc cos 4 9. 

momentum flux than in models A and B due to the high col- 
limation. The TS is also confined more tightly by the ISM 
ram pressure. 

3.4 Split-monopole versus wave-like wind 

In principle, the results above afford a means of distinguish- 
ing between wind models with different momentum flux dis- 
tributions. In reality, however, the conclusion from section 



Figure 8. y-z section of model Al. The colors denote the den- 
sity in logarithmic scale, as defined in figure El The symbol m 
the top left indicates the projected wind symmetry axis and 
the momentum flux (polar plot), as in the caption of figure 121 
The solid curve shows the theoretical solution <A2I . scaled to 
match the CD (dashed curve). This model describes a wind whose 
momentum flux is concentrated in the equatorial plane, with 
g w (6) = sin 2 (9 + 7/3 



13.31 and figure [5] that the shapes of the CD and BS de- 
pend weakly on g w (9), seriously hampers this sort of ex- 
pe riment. For example , the w ave-l ike dipole wind analyzed 
by iMelatos fc Melrosel Jl99fJ) and iMelatol lll997l) predicts 
g w iff) = cos 2 8 + 1 (for a point dipole, as calculated in 
appen dix ITU , while the split-monopole wind analyzed by 
iBogovalo ^19991) predicts g w {0) = sin 2 + 1/cr, yet these 
models produce very similar shapes for the CD and BS, as 
we now show. 

Consider an equatorially dominated wind, C2 = —1, im- 
plying g w {9) = 7/3 + sin 2 8. This sort of wind is favo red by 
some PWN observations ijChatteriee fc Cordesl EoO^. Fig- 
ure |H| shows the y-z section of model Al. We can hardly see 
any characteristic features in the bow shock. The CD is still 
described by the analytic solution. The TS looks similar to 
the TS in model A, mirrored along the x-z plane. Indeed, 
this simulation strongly resembles the situation of a pole- 
dominated flux with A = 135° (cf. A = 45° in model A). 
We conclude that it will be hard to distinguish these cases 
observationally. 

Out of interest, a comparison between the y-z section 
and x-z section of these models can, in principle, reveal the 
underlying wind structure. Such a comparison is depicted 
in figure |U] The pole-dominated wind (left panel) has an 
anisotropic TS, while the TS in the equatorially-dominated 
wind (right panel) resembles a completely isotropic out- 
flow. Unfortunately, such a comparison requires orthogonal 
lines of sight and is not accessible by observation. The com- 
putation of synchrotron emission maps for both situations 
from different angles may provide additional information, 
if Doppler boosting and column density effects break the 
degeneracy in figure [S] However, reliable synchrotron maps 




Figure 9. x-z section of model A (left) and model Al (right). The colors denote the density in logarithmic scale, as defined in figure 01 
The symbol in the top left indicates the wind symmetry axis, as in the caption of figure|5] The solid curve shows the theoretical solution 
IA2I . scaled to match the CD (dashed curve). 



require RMHD simulations and lie outside the scope of this 
paper. 



4 COOLING 

4.1 Numerical implementation 

The analytic solution derived bv lWilkinl Jl99fil200o|) is exact 
in the single-layer, thin-shock limit, where both the shocked 
ISM and the wind cool radiatively (via spectral line emis- 
sion and synchrotron radiation respectively) faster than the 
characteristic flow time-scale. In this section, we quantify 
how rapidly such cooling must occur for the thin-shock ap- 
proximation to be valid. We also present some simulations 
of shocks where cooling is not efficient and draw attention to 
h ow the structur e of these shocks differs from the predictions 
of lWilkinlfeOOOT) . 

Flash implements optically thin cooling by adding a 
sink term — A(x,t) to the right-hand side of equation J3J|. 
Radiation is emitted from an optically thin plasma at a 
rate per unit volume A(T) = Xo{l — Xo)p 2 rrip 2 P(T), where 
p is the density, Xq is the ionisation fraction, m p is the 
proton mass, and P( T) is a loss function that depends on 
the temperature only ([Raymond et alJll976l : ICox fc Tuckerl 
119691: iRosner et alJll97Sl) . In our temperature range, 10 4 K 
^ T ^ 10 5 K, hydrogen line cooling dominates and we can 
write P(T) ~ 10 -22 erg s _1 cm 3 . The postshock tempera- 
ture Tps can be extracted directly from the simulation or 
estimat ed from the Bernoulli equatio n and entropy conser- 
vation feucciantini fc Bandierall200ll) . With the latter ap- 
proach, the result is feTps = 3m p u 2 /32, where ks is Boltz- 
mann's constant. 

For a slowly moving pulsar, like PSR J2124— 3358 which 
has v p = 61 km s" 1 and p m = 1.6 x 10~ 24 g cm -3 



JCaensler et al.ll2002ft . we find T PS = 4.2 x 10 4 K and 

A = 6.1 x 10" 21 erg s" 1 cm" 3 , if we assume a moderate 
ionisation fraction Xo = 0.2. The internal energy per unit 
volume in the postshock flow is e — ( r y—l)~ 1 (N a kB/A)Tp = 
1.4 x 10" 11 erg cm -3 where A — 1 g is the average molec- 
ular mass of the hydrogen plasma, from which we deduce 
a cooling time r = e p /A = 2.3 x 10 9 s. In this object, 
the optical bow shock extends 65"or / = 2.3 x 10 18 cm, 
implying a flow time tf = l/v p = 3.8 x 10 11 s and hence 
r <C tf. This example shows that cooling is important not 
only in high-luminosity pulsar s and in a high density ISM 
iBucciantini fc Band icra 200 lj), but also if the shock is suf- 
ficiently extended, i.e. if the typical flow time, depending on 
v p , approaches the t ypical cooling time, which is a function 
of v p and p m . While lBucciantini fc Bandieral i200ll) argued 
that the typical cooling length scale is the width of the BS 
plus the distance between BS and CD, we include the whole 
visible BS, allowing a longer time for the shocked ISM to 
cool. This is justified because the material in the post-BS 
region continues to cool and participate in the dynamics of 
the system. 

We can perform a similar calculation for a fast pulsar, 
like PSR J1747-2958 (the Mouse) w hich has v v = 600 km 
s" 1 and p m = 5.0 x 10" 25 g cm" 3 dCaensler et alJl2002T) . 
corresponding to T PS = 4.1 x 10 6 K and A = 1.7 x 10" 24 erg 
s" 1 cm" 3 . The cooling efficiency here is lower than in the 
previous example, since hydrogen line co oling is slower at 
high temperatures iRavmond et al.|ll97fj) and p m is lower. 
The cooling time r = 1.5 x 10 14 s is long compared to the 
flow time-scale tf = 1.5 x 10 10 s; the length of the BS is 
/ = 0.3 pc. 

To treat Ha cooling properly, one must solve the Boltz- 
mann equation for the neutral atom distribution as a func- 
tion of position, including the ion-neutral reactions, a cal- 
culation which lies beyond the scope of this article. Instead, 
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in order to quantify the overall effects of cooling and ver- 
ify the validity of the analytic thin-shock approximation, we 
parameterize the cooling function as 

A(x,t)=Ao/pe(x,t), (9) 

with Ao = 10~ 7 s _1 , such that the dimensionless parame- 
ter / controls the local cooling time through r = (Ao/) -1 
(which is then spatially uniform by construction, a crude 
approximation) . 

4.2 Shock widths 

Figure [T01 shows the y-z sections of models CI to C4. These 
four models are identical (g w oc cos 2 9, A = 35°; see ta- 
ble □ , except that / varies from 1CT 4 (CI) to 10° (C4). 
For more efficient cooling, the BS becomes thinner (compare 
upper- left and bottom-right panels in the figure). Similarly, 
the TS becomes more extended with increasing cooling ef- 
ficiency, until the wind cavity fills up the BS interior and 
the TS becomes indistinguishable from the BS for highly 
efficient cooling (bottom-right panel). Since the pressure in 
the PS region approximately equals the thermal ISM pres- 
sure, and since efficient cooling reduces the PS pressure, the 
wind ram pressure drives the termination shock outward as 
/ increases, as seen in figure 1101 On the other hand, the 
thickness of the bow shock is controlled by the ratio be- 
tween the thermal and ram pressures in the ISM. The BS 
region therefore becomes thinner for more efficient cooling, 
as we pass from model CI to C4. For / = 1, the TS and BS 
almost touch and the one-layer, thin-shock analytic solution 
given by equation <A21 applies as in the lower-right panel of 

figure run 



5 SMOOTH ISM DENSITY GRADIENT 

Density gradients in the ISM have been invoked to 
explain the asymmetric features observed in the bow 
shocks around PSR B0740-28, PSR B 2224+65 (the Gui- 
tar nebula), and PSR J2124-3358 llJones et alj 12001 
IChatteriee fc CordeJ2002t I Gaensler et al J|200j). The PWN 
around PSR B0740— 28 is shaped like a key hole, with a 
nearly circular head and a divergent tail, starting ~ 20i?o 
downstream from the apex. The Guitar nebula consists of a 
bright head, an elongated neck, and a limb-brightened body. 
The PWN around PSR J2124-3358 exhibits an asymmetric 
head and a distinct kink ~ 10i?o downstream from the apex. 

5.1 Numerical implementation 

To explore the effect of a smooth density gradient, we carry 
out a set of four simulations (D-G in table in which the 
ISM has an exponential density profile 

Pm(y) = p m ,o exp(-y/H), (10) 

where H is a length scale, and p m ,o substitutes for p m in 
table Q We generally choose H ~ Rq ~ 0.01 pc, consistent 
with the wavelength o f turbulent fluctuations in the warm 
ISM JPeshpandel2000l) . The gradient is perpendicular to the 
pulsar's motion. The density profile is initialized and then 
maintained at later times by setting p m (y) at the inflow 
boundary at z = z max to p m (y) = p m ,o exp(— y/Ro), and 



choosing p m and e m consistently. The time-scale over which 
the associated pressure gradient smoothes out p, H/c m ~ 
10 10 s, exceeds the flow time-scale, tf = (z mK — z m i„)/D ra ~ 
6 x 10 9 s, so the density gradient remains nearly constant. 

The analytic solution for the shape of a thin bow shock 
in an exponential density gradient is given by equation 
JAT21 . 



5.2 Tilted shock 

Figure 1111 shows the y-z sections of models D-G, together 
with the analytic solution for the BS (solid curve) and the 
CD (dashed curve). In these models, the exponential den- 
sity profile changes from shallow (H — L5.Ro for model D, 
upper-left panel) to steep (H = 0.25-Ro for model G, bottom- 
right panel). One distinctive feature is that the low-density 
side (y > 0) of the shock structure, where the CD, TS, and 
PS are situated, is tilted with respect to v m by an angle of 
up to 180° for H/Rq = 0.25. This occurs because the ISM 
ram pressure, p m u^, is dominated by the wind ram pressure 
for y > 0, pushing the CD further away from the pulsar. In 
fact, this part of the CD is appro ximated well by the ana- 
lytic formula SKHi (IWirkinll200ol) . The BS separates from 
the CD and broadens as the thermal pressure of the ISM 
in the low-density region increases relative to the ISM ram 
pressure. 

On the high-density side (y < 0) of the bow shock, the 
opposite situation prevails. Here, the ram pressure of the 
ISM exceeds its thermal pressure, so that the BS is well 
approximated by the thin-shell solution. The CD curves to- 
wards the lower density region, pushed by the thermal ISM 
pressure. The opening angle of the CD increases as H de- 
creases: we find 40°, 60°, 103°, and 124° for models D, E, 
F, and G respectively, as the ISM ram pressure (for y > 0) 
decreases with increasing y. 

Applying the above results to observations of PWN, 
we expect the Ha surface brightness (oc column density) to 
mimic the volume density contours as in Figure 1111 espe- 
cially where the thickness of the BS increases from y < to 
y > 0. Such a variation in the Ha flux has not been observed 
in the PWN aroun d PSR B0740-28 PSR B2224+65, and 
PSR J2124— 3358 lOones et alj Eool IChatteriee fc Cordesl 
l2002UGaensler et al]|2002ft . It therefore seems unlikely that 
a smooth density gradient is responsible for the peculiar 
morp hologies observed in these objects llChatteriee et all 



6 WALL IN THE ISM 

The IS M is inhomogenou s on length scales from kpc down 
to AU JPeshpandeil2000D . Hence a fast pulsar is likely to 
encounter singular obstacles or barr iers, such as the edges 
of O-star bubbles dNaze et al.l 12002?) . which effectively act 
as "walls" . A wall was invoked to explain t he kink observed 
in th e nebula around PSR J2124— 3358 llChatteriee et all 
l200rJ) . U nlike when a pulsar interac ts with its supernova 
remnant llvan der Swaluw et al.ll2003T) . the pulsar does not 
always strike the wall head on, giving rise to a highly asym- 
metric interaction. 
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Figure 10. y-z section of models CI (/ = 1CT 4 , upper left), C2 (/ = 10" 3 , upper right), C3 (/ = 10~ 2 , bottom left), and C4 (/ = 1, 
bottom right). The colors denote the density in logarithmic scale as indicated by the color bar (g cm" 3 ), while the solid curve shows 
the analytic solution IA2I . The symbol in the top left indicates the projected wind symmetry axis and momentum flux (polar plot), as 
in the caption of figure |2l For more eff icient cooling, e.g. C4, the TS lengthens and the BS becomes thinner, until the TS and BS are 
no longer separated. The IWilkirJ fcOOfj) solution (solid curve) matches the bow shock nearly exactly (bottom right), with a deviation of 
< 8% from the actual BS position. Here, the bumps at the BS are an artifact of the discrete grid. They are not seen in C1-C3, where 
small density fluctuations are smoothed by the thermal pressure. However, C1-C3 exhibit strong KH instabilites at the BS. 



6.1 Numerical implementation 

In order to model the wall, we introduce the (signed) co- 
ordinate w along the direction defined by cos ay — sinaz, 
perpendicular to the wall plane. In model Wl, the wall is 
infinitely extended (i.e., a ramp). Its profile can be written 



In models SI and S2, we truncate the wall on one side. Its 
profile can be written as 

p(w) = p m + (pwaii - p m )[l + exp(a7/A r )]~ 

+0(— w - 8A r )(p wa n - p m ) exp[(uj + 8A r )/A c j], 

(12) 



p{w) = (pwall - Pm) [1 + exp(lI7/A r )] 1 + p„ 



where 6(zo) is the usual Heaviside function, denned as 
(11) 6(zu) = {0(vj < 0), l(ro ^ 0)}. In this notation, w = 
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Figure 11. sections of model D (H = 1.5-Ro, upper left), E (H = l.ORo, upper right), F (H = 0.5R , lower left), and G (H = 0.25R , 
lower right). The colors denote the density in log 10 (g cm -3 ); note that the color scale differs from that in the other figures. The solid 
curve shows the theoretical solution IA12I . scaled to match the CD (dashed curve). The analytic approximation clearly breaks down for 
models D— G (see text). The BS deviates from the CD where the ISM ram pressure exceeds the thermal pressure. Axes are labelled in 
units of R = 1.85 x 10 16 cm. 



denotes the center of the smooth rise in equations HI IB and 
ill 2t . The maximum density is p wa ,u = 1.30 x 1CP 3 g cm -3 . 



6.2 Kinks in the bow shock 

In this section, we examine the effect of a wall on the shock 
structure and work out how PWN observations probe inho- 
mogenities in the ISM. We consider both a ramp-like [equa- 
tion pip] and a truncated [equation 112^ ] wall. 

Figure 1121 shows a simulation (model Wl) where the 
pulsar hits a ramp-like wall with a = 26° and density con- 



trast rj = pwaii I pm = 2.5. The wind is asymmetric, with 
A = 90°and ft tilted towards the observer {cj> = 45°), result- 
ing in a slightly asymmetric PS flow. The top and bottom 
rows of figure 1121 depict snapshots of x-z and y-z sections 
respectively. The y-z sections demonstrate the asymmetry 
best. On impact, the BS is compressed by the increased ram 
pressure of the wall. A kink appears at the transition point 
between the wall and the ambient ISM. Both the gradient 
in Ha brightness and the kink should be observable in prin- 
ciple. As the pulsar proceeds into the high-density region 
(see middle panels of Figure 1120 , the CD approaches the 
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Figure 12. Time series of x-z (top row) and y-z (bottom row) sections of model Wl. The columns are at times t = l.7tf (left), t = 2.5tf 
(middle), and t = 3Atf (right), where tf = 6 X 10 9 s is the time for the ISM to cross the simulation volume. The colors denote the density 
in log 10 (g cm -3 ), while the symbol in the top left indicates the projected wind symmetry axis and momentum flux (polar plot), as in the 
caption of figurc[2] The density of the wall is 4 times higher than the ISM density. The momentum flux is anisotropic (A = 90°) and f2 is 
tilted towards the observer (0 = 45°). The bow shock is not disrupted by the wall and shows characteristic kinks and Kclvin-Hclmholtz 
instabilities. 



pulsar due to the increased ISM ram pressure. Likewise, the 
ISM thermal pressure increases with increasing p m and, by 
the argument given in section 13.31 the rear surface of the 
TS comes closer to the pulsar. Just as for a smooth ISM 
density gradient (cf. Figure II It , the shocks are tilted to- 
wards the low density side. When the pulsar proceeds into 
the densest part of the ISM, the BS, CD, and TS approach 
even closer to the pulsar until they finally reach the state in 
the rightmost panels. 

In this stationary state, KH instabilites triggered by 
perturbations during the collision with the wall occur at the 
right hand surfaces of the CD. In the absence of gravity, the 
KH instability occurs for all wave numbers and hence on all 
length-scales down to the grid resolution (~ 0.1i?o in this 
simulation). However, the collision with the wall introduces 
a velocity perturbation that is large (~ v m ) compared to 
the numerical noise. This explains why we see the instability 
in the right panels but not the left panels, where the wall 
has not yet reached the BS. Indeed, the length scale of the 
density rise due to the wall agrees with the wavelength of 
the observed KH instability (~ Ro). 

In the lower middle panel, the CD downstream from the 
apex is slightly lopsided towards the higher density (left) 
side. This occurs because the PS flow from the side sur- 
faces of the TS travels faster than the ISM. The region in- 



side the CD therefore loses pressure support by the shocked 
wind material, and the shocked ISM material pushes into 
the space. The same effect can be seen in the lower right 
panel on the lower density (right) side, where the KH insta- 
bilities favor the right-hand side. The top row {x-z sections) 
basically resembles a head-on collision, similar to the in- 
teraction of a PWN with th e pulsar's supernova remnant 
llvan der Swaluw et al.lEo03l) . Note also the broadening of 
the lower part of the BS (top right panel), which demon- 
strates the delay (~ 5Ro/v m ) before the shock adjusts to 
the higher external density p m . 

The Ha emission eman ates chiefly from the region be - 
tween the BS and the CD llBucciantini fc Bandieral 120011) . 
which is only moderately affected by the passing wall. Fur- 
thermor e, ISM inhomogenities often consist of highly ionised 
matter JChatteriee et alj|200rj) . so that the Ha luminosity 
observed as the wall passes is dominated by the afterglow of 
the neutral part of the non-wall ISM. Negligible brightness 
variations are observable in either case. Nevertheless, as we 
discuss in section |7| a comprehensive treatment of the emis- 
sion requires the solution of the Boltzmann equation for the 
neutral and ionized species, a task beyond the scope of this 
article. 

In models SA (a = 26.6°) and SB (a = 63.4°), dis- 
played in figure tTSI the density ramp is replaced by a wall of 
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Figure 13. x-z (top row) and y-z (bottom row) sections of model SA (left column) and model SB (right column). The colors denote the 
density in log 10 (g cm -3 ), while the symbol in the top left indicates the projected wind symmetry axis and momentum flux (polar plot), 
as in the caption of figure |2 In these simulations, the pulsar travels into a wall with a finite width. The momentum flux is anisotropic 
(A = 90°) and f2 is tilted away from the observer (0 = —45°). 



finite width, which is truncated as in equation 112H . The top 
and bottom row again shows x-z and y-z sections respec- 
tively. As above, the x-z sections resemble a head-on colli- 
sion. The base of the BS broadens in the low-density region 
and narrows at the density maximum. Above the wall, the 
BS readjusts to the lower ISM density. Due to the smooth 
gradient (Ad S> A r ), the BS looks "egg-shaped", an effect 
that can been seen more clearly in the upper-right panel, 
where the inclination of the wall is lower (a = 64°). The 
kink is also more distinct, as seen in the y-z sections (lower 
panels), because the BS expands in the low density region 
above the wall. The CD and TS approach the pulsar for the 



same reasons discussed above, but the smooth gradient is 
responsible for the bulbous shape. Again, as in figure [T21 we 
can see in the lower-right corners of the lower panels how the 
shocked wind material inside the CD loses pressure support 
and the shocked ISM material diffuses in. 



7 DISCUSSION 

In this paper, we perform three-dimensional simulations of 
pulsar bow shocks in the ISM under a variety of asymmetry- 
inducing conditions: an anisotropic pulsar wind, an ISM den- 
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sity grad ient, and an ISM wal l. We validate the analytic so- 
lution of IWilkinl lll996l . |200Cf) for an anisotropic wind and 
show that it remains a reasonable approximation for the 
thick-shock case, with a deviation of ~ 10% within ~ 5Ro 
of the apex. We show that the basic multilayer shock s truc- 
ture (TS, CD and BS) o btained bv IBucciantini <|2002F l and 
iBucciantini et all (120051) for an isotropic wind is preserved 
in an anisotropic wind as well. We find that the thermal 
pressure smoothes the BS. The underlying anisotropy of the 
momentum flux thus remains concealed, making it difficult 
to gauge the ISM density from the stagnation point, since 
generally neither the distribution of the momentum flux nor 
the orientation of the pulsar spin axis is known. 

Contrary to previous estimates 

Bucciantini fc Bandieral 1200 J) , we show that cooling 
can be important for slow pulsars in a low-density ISM, 
such as PSR J2124-3358. This is because the shocked ISM 
continues to take part in the dynamics of the system and 
thus influences the overall shape. There are observational 
consequences: cooling tends to increase the overall size of 
the TS while the BS becomes thinner. Naively, a thinner 
BS is expected to be brighter due to limb brightening. 
However, in order to quantify the effects of cooling on the 
optical emission, it is necessary to solve the collisionless 
Boltzmann equation for the neutral hydrogen atoms, a task 
beyond the scope of this paper. 

Our simulations are non-relat ivistic and do not i nclude 
the effects of a magnetic field. IBucciantini et all (120051) 
showed that a proper relativistic treatment increases the 
KH instabilities arising at the CD by enhancing the velocity 
shear at the interface. In principle, therefore, shocked ISM 
material can contaminate the wind material, altering the in- 
ner flow structure and perhaps modulating the synchrotron 
flux in time, although no simulation has demonstrated the 
latter effect to date. A magnetic field affects the flow more 
drastically, in two ways. First, the standoff distance increases 
because the magnetic field adds to the pressure at the lead- 
ing surface. This i s a ru naway process in the simulations of 
IBucciantini et all l)2005|) . who neglected resistivity. Second, 
the TS becomes convex because the magnetic pressure is 
not uniform. The azimuthal magnetic field builds up near 
the pulsar, increasing the magnetic pressure. This effect is 
strongest at z v and decreases as z decreases. 

An interesting question is whether the magnetic field 
in PWN can produce the torus-ring str ucture found by 
iKomissarov fc Lvubarskvl l)200.^ . l20o"iF and lDel Zanna et all 
|120mF ~ where a jet is formed by circulating back flow of 
the shocked wind material colliminated subsequently by the 
magnetic field. This circulation takes place outside the TS 
and is suppressed by strong confinement in PWN. How- 
ever, a detailed understanding of the magnetic field structure 
when the momentum flux is asymmetric must await future 
simulations. 

PWN are commonly detected as radio or X-ray syn- 
chrotron sources. Charged wind leptons are accelerated at 
the termination shock and subsequently gyrate in the local 
magnetic field, emitting synchrotron radiation. Typically, 
the emission is observed from three different zones: the TS it- 
self, the PS region, and the head of t he BS. For an identifica - 
tion of these zones in the Mouse, see lOaensler et all l)2004l . 
A further theoretical study of the synchrotron emission will 
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be u ndertaken in a forthcoming paper iChatteriee et alJ 
120061) . Here, we simply make a few qualitative points. 

(i) Generally, we observe a variety of TS shapes from in- 
clined tori (Crab Nebula, where the emission is equatorially 
concentrated) to cylinders (as seen in the Mouse, for a nearly 
isotropic wind) . The emission from the TS depends strongly 
on the anisotropy of the wind. The more peculiar the shape, 
the more information it can provide about the inclination of 
the spin axis with respect to the direction of motion as well 
as the angular distribution of wind momentum flux. The 
asymmetry in the TS (figures |5] and directly corresponds 
to the axis of maximum emission. 

(ii) For isotropic wind emission, one expects the PS flow 
to be observed as an uncollimated X-ray tail feucciantinil 
I2002T) . However, if the wind is anisotropic (cf. Figures 
Q, our simulations predict overluminous and underluminous 
regions (corresponding to high and low densities) in the PS 
flow, even two separated tails in extreme cases (as in figure 

(iii) The observed appearance of the BS is strongly in- 
fluenced by the emission from particles accelerated at the 
head of the BS and subsequently swept back along the 
CD. Since only the flow component normal to the shock 
slows down, the flow in the swept-back region is generally 
trans- relativistic (~ c). If this part of the PS flow [region 
Bl in IChatteriee et all teOOfJl l is laminar along the CD, 
Dop pler beaming may rende r the emitted radiation invis- 
ible (IBucciantini et alll2005l) . An anisotropic wind, how- 
ever, may have a sufficiently large velocit y component to- 
ward s the line of sight to allow detection ijChatteriee et all 
120061) . In the case of a jet-like wind, the BS can exhibit a 
one-sided X-ray tail whose receding half is suppressed by 
Doppler beaming as in the PWN around PSR J2124-3358 
dChatteriee et alj|200rj. whose model B corresponds to our 
model SA). 

Some PWN are also detected in neutral hydrogen emis- 
sion lines JCaensler fc Slanell200ef) . Neutral hydrogen atoms 
are excited collisionally or by charge exchange and then de- 
excite radiatively in the region between the BS and CD. 
Hence, the Ha luminosity depends crucially on v p and Xo- 
A proper kinetic (Boltzmann) trea tment incorporating ion- 
neutra l reactions is discussed by IBucciantini fc Bandieral 
feOOlf) . It is clear that purely hydrodynamical simulations 
like ours cannot predict the observed surf ace brightness re- 
liably . However, in a paper in preparation llChatteriee et all 
1200a) . we show that the overall shape of the bow shock is 
predicted reliably. 

If high-resolution observations of the shock apex can 
resolve the characteristic kink in the BS at the latitude 
where the momentum flux is a minimum, they will shed light 
on the angular distribution of the pulsar wind momentum 
flux and hence the electrodynamics of the pulsar magne- 
tosphere ormssarov & Lvubarskvl 1200.1 

12004 iMelatos fc Melrose! Il99(t lMelatoslll997l . 120021) . How- 
ever, it is generally difficult to distinguish between equatori- 
ally dominated (e.g., split-monopole) and pole-dominated 
(e.g., wave-like dipole) wind models, beacuse the shock 
structure is degenerate with respect to several parameters 
and does not depend sensitively on the momentum flux dis- 
tribution (see sections 13.31 and 13.40 . Only for favored incli- 
nations, e.g. when fl points directly towards us and the TS 
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ring is resolved (as in the Crab nebula), can the synchrotron 
emission from the TS be used to make more definite state- 
ments. Again, RMHD simulations for anisotropic winds and 
synchrotron emission maps are needed to break the degen- 
eracies. 

Pulsar bow shocks probe the local small-scale structure 
of the ISM. An ISM density gradient modi fies the shape 
of the bow shock substantially away from the lWilkinl l|200(t) 
solution. Since Ha emission is generated behind the BS, den- 
sity fluctuations with a length-scale > Ro should be clearly 
visible. If the gradient is nearly perpendicular to v m , we ex- 
pect the surface brightness of the nebula to be asymmetric 
as well, even one-sided in extreme cases. In addition, the 
TS is tilted towards the low density side, producing asym- 
metric synchrotron emission. Steep density gradients with 
length-scale < Ro (i.e. walls) produce characteristic kink- 
like features. For example, a wall consisting of mainly ionized 
matter can account for the obs erved shape of the PW N asso- 
ciated with PSR J2124-3358 ijChatteriee et alJl2006h . If the 
pulsar encounters a truncated wall with a finite width, the 
bow shock develops a peculiar shape: an egg-like head and 
a divergent tail (as in figure 1131 . This shape is strikingly 
similar to the Ha emission observed from PSR B0740— 28 
J Jones e t al.||2002| ). Multi-epoch ob servations of the Guitar 
nebula JChatteriee fc Cordesl 120041) reveal time-dependent 
behaviour of this PWN that corresponds to the overall situ- 
ation shown in figure (although not to the time between 
the snapshots in figure 1121 of course, which are separated 
~ 10 2 yr). When the pulsar runs into the higher density 
region, the wind is confined more strongly and the BS is 
pushed closer to the pulsar, explaining the narrow head and 
diverging tail. Furthermore, we note that the Guitar neb- 
ula is significantly brighter near the strongly confined head. 
This can be explained in terms of the increased density of 
the shocked ISM, if we assume that the high density region 
carries a neutral fraction comparable to the low density re- 
gion. On the other hand, our simulations do not exhibit the 
rear shock seen in the Guitar nebula, which could be due to 
a higher r\ than considered in models SA and SB. 

Summarizing our results, we conclude that the 
anisotropy of the wind momentum flux alone cannot ex- 
plain odd bow shock morphologies. Instead, it is necessary 
to take into account external effects, like ISM density gra- 
dients or walls. Conversely, because the shape of the bow 
shock is degenerate with respect to several parameters, it is 
difficult to infer the angular distribution of the wind momen- 
tum flux from the Ha radiation observed. Highly resolved 
radio and X-ray observations, combined with synchrotron 
emission maps from RMHD simulations, may break these 
degeneracies in the future. Finally, we caution the reader 
that the results presented here explore a very limited por- 
tion of the (large) parameter space of the problem. 
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APPENDIX A: ANALYTIC BOW SHOCK 
FORMULAS IN THE THIN-SHELL LIMIT 



l 9' + w cos 9' . 



(A3) 



IWilkinl dl996T) showed that the shape of the CD when the 
pulsar wind momentum flux is isotropic can be written as: 

R(9) = R o csc0[3(l - 6>cot6>)] 1/2 (Al) 

in the limit of a thin shock. Here, R(9) gives the distance 
from the pulsar position to the bowshock, 9 denotes the 
colatitude measured relative to the orientation of the sym- 
metry axis (see section |3J, and Ro is the standoff distance 
defined in equation JSJ . 

The formalism was extended to an anisotropic pulsar 
wind with p = 2, by I Wilkin] feOOOT) . who found 

R(9) = i? csc6»'{3(l-6i'cot6l')[co + ^(3D 2 +'u; 2 )] 

+ ^j-[(v 2 - w 2 ) sin 2 9' + vw(29' - sin 26»')]} 1/2 , 

(A2) 

with v — sin 0' sin A, w = cos A, and Co = 1 — C2/3, if the 
momentum flux is of the form p w {9) = po, w (co + C2 cos p 9) as 
in equation (|7J . 9' and <j>' are the usual spherical coordinates 
with respect to the z axis, such that 



Note that this notation differs from IWilkinl (l200fj) : his (9* 
corresponds to our 9, while his 9 and (j> correspond to our 9' 
and (/>' respectively. 

We derive a formula for R(9 ) in th e case p — 4 by 
applying the formalism of IWilki nl j200C j). The case p = 4 
is relevant to a collimated jet JChatteriee et alJl200rJ) or an 
extended vacuum dipole (see appendix [HJ. In side the wind 
cavity, the momentum flux can be written as JWilkinll2000Tl 



p w (r,9)vi = E 
2nr z v v 

with 

gw{9) = C0 + C2 cos 4 



(A4) 



(A5) 



Here, the overall energy loss E is related to the simulation 
parameters po,u> and r w through E — 2Trr 2 ^po, w v^ 1 . The nor- 
malisation 2 



d(j} / d9 sin 9 g w {9) = 4tt 



(A6) 



requires Co = 1 — C2/5 for p = 4. The incident momentum 
flux from the wind on the shell can then be written as 



') = 



E 



■Gw(9 



(A7) 



Here, G w — G w ^zb + G WtZ z, where the cylindrical coor- 
dinates zu and z are denned with respect to the symmetry 
axis, is a dimensionless function given by 



G w = I d9 g w (zb sin 9 + z cos 9 ) sin 



(A8) 



From ifAlJ and ljA"5|l . we can compute g w and hence G w as 
a function of 9' and <j>'\ 

g w = co + C2(sin9'v + cos9'w) 4 . (A9) 

It can now be shown JWilkinl2000D that the bow shock shape 
is described by 

R(9',4>') = R csc9'[-6(G w ^cot9' - G w , z )] 1/2 . (A10) 
Substituting 1A8H and IA9I into IAIOII finally yields 

R = — esc 9'{- csc 9'{c [1929' cos 0' - 192 sin 9'] (All) 

8 

+c 2 [«™ 3 (16cos6»' - 12 cos 36' -4 cos 56' -966' sin 9') 
+v 3 w{32cos9' -36 cos 36>' +4cos56>' -966' sin 9') 
+v 2 w 2 (14461' cos 6' - 168 sin 9' + 18 sin 36' - 6 sin 56') 
+v 4 (1209' cos 9' - 80 sin 0' - 15sin36»' + sin 59 ) 
+w 4 (249' cos 9' - 56 sin 6' + 9 sin 36' + sin 56»')]}} 1/2 . 

Figure |A"T1 shows the two solutions 1A2I and HA 1 111 for 

A = 7r/4, co = 0, and C2 = 3. Although the shapes are nearly 
identical, the p = 4 surface is slimmer. The difference is seen 
more clearly in figure lA"2l which shows the y-z section of the 
bow shock. The p = 4 case (dashed curve) touches the p = 2 

2 Although this normalisation is not necessary for the analytic 
solution, it is required for equation J7J to hold. If it is not fulfilled, 
the definition of the scale factor Ro changes accordingly. 
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Figure Al. Comparison of the bow shock shape for p = 2 (left) 
and p = 4 (right), given A = 7r/4, co = and C2 = 3. We note 
the strong similarity between these solutions, the only difference 
being the slimmer shape for p = 4. 
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Figure A2. y-z section of the bow shock for p = 2 (solid curve) 
and p = 4 (dashed curve). The strong similarity hampers efforts 
to distinguish between the two cases observationally. 

shock (solid curve) at the latitude of maximum momentum 
flux, 9 = 7r/4 (upper-right corner). Both surfaces show a 
characteristic kink at the latitude where the momentum flux 
is a minimum. 

By a similar calculation, we can derive the shape of the 
BS for an exp onential density gradient perpendicular to v m 
is 



Ha tor an exp onential den 
JWilkinll2000l) . The result 

R(8, cj>) = H sec 4> sin -1 



ay, Y), 

where y s is the solution of 

- [y - 2 + exp(-y)(y + 2)] = JyL 
y 6 

and t/ax is the normalised standoff distance, 

Z/ax^^ COS ^[3(1-0 COt 0)] 1/2 . 
12 



(A12) 



(A13) 



(A14) 



APPENDIX B: ANGULAR DISTRIBUTION OF 
THE MOMENTUM FLUX FOR AN EXTENDED 
VACUUM DIPOLE 

Although the electrodynamics of the pulsar magnetosphere 
and wind remain unsolved, it is commonly assumed that the 
angular distribution of the energy flux at the light cylinder 
is preserved in the transition to a kinetic -energy-dominated 
outflow at the TS (SeT Zanna et al .120041) . In this appendix, 
we calculate the far-field Poynting flux for point-like and 
extended dipoles rotating in vacuo. 

The relevant electromagnetic far field compo nents are 
given by the real parts of the following expressions iMelatosI 



B 



B 



[a(xo) — b(xo)] sin a cos 8e +T> + 0(x 2 ) 

[a(x ) cos 29 - b(x )] sin ae i( * +v) + 0(x~ 
-B e ). 



, (Bl) 

2 ), (B2) 
(B3) 



ijh 
2x 
_Bo 
2x 

(Eg, E$) = c(B^, 

In JB 11 1— dT33ll . x = cor/c denotes a normalized radial coor- 
dinate, uj is the pulsar's angular velocity, xo — lotq/c is the 
normalized effective radius of the corotating magnetosphere, 
rj = 4> — wt is the phase of the electromagnetic wave, and Bo 
is the magnitude of the stellar magnetic field at the poles. 
Constants a and b are given by 



J2 
Xq 



x h' 2 (x ) + h 2 (x ) 



and 
b = 



.'•<> 



(B4) 



(B5) 



hi(x ) 

where hi (x) and /12 (x) are the first and second-order spher- 
ical Hankel functions of the first kind. The r component of 
the Poynting flux can then be computed from 

1 

2/LtO ' 

= So(\a — b\ 2 cos 2 9 + \a cos 29 ■ 

with So = (c/2^ )(Bq/2x) 2 . 

We can simplify 1B6II using the explicit form of the Han- 
kel functions, 



(E x B*) 



6| 2 ) sin 2 a 



(B6) 



hi(x) =e lx (-- - \ 
X x z 



and 
h 2 {x) 



3 ./i 
\x 



(B7) 



(B8) 



In the case of a point dipole, we have xq 1. Retaining the 
leading order terms, we obtain 



1 5 — ixn 

6 

. 3 — ixQ 

-ix e 



and hence 

S r = So2;o(cos 2 9 + 1) sin 2 a, 

which reduces to the case p — 2, 
fcOOOl). 



Co 



(B9) 
(BIO) 

(Bll) 

= C2 = 1 in IWilkinl 



We also calculate the Poynting flux for an extended 
dipole, as proposed by iMelatosl (Il997l) . For the special case 



xo 



1, an explicit calculation yields 
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(2cos 4 + 6cos 2 + 5)sin 2 a. (B12) 



